home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Fritz: All Fritz
/
All Fritz.zip
/
All Fritz
/
FILES
/
PROGMISC
/
PCSSP.LZH
/
PC-SSP.ZIP
/
MATINSL.ZIP
/
DMLSS.FOR
< prev
next >
Wrap
Text File
|
1985-11-29
|
6KB
|
203 lines
C
C ..................................................................
C
C SUBROUTINE DMLSS
C
C PURPOSE
C SUBROUTINE DMLSS IS THE SECOND STEP IN THE PROCEDURE FOR
C CALCULATING THE LEAST SQUARES SOLUTION OF MINIMAL LENGTH
C OF A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH SYMMETRIC
C POSITIVE SEMI-DEFINITE COEFFICIENT MATRIX.
C
C USAGE
C CALL DMLSS(A,N,IRANK,TRAC,INC,RHS,IER)
C
C DESCRIPTION OF PARAMETERS
C A - COEFFICIENT MATRIX IN FACTORED FORM AS GENERATED
C BY SUBROUTINE MFSS FROM INITIALLY GIVEN SYMMETRIC
C COEFFICIENT MATRIX A STORED IN N*(N+1)/2 LOCATIONS
C A REMAINS UNCHANGED
C A MUST BE OF DOUBLE PRECISION
C N - DIMENSION OF COEFFICIENT MATRIX
C IRANK - RANK OF COEFFICIENT MATRIX, CALCULATED BY MEANS OF
C SUBROUTINE DMFSS
C TRAC - VECTOR OF DIMENSION N CONTAINING THE
C SUBSCRIPTS OF PIVOT ROWS AND COLUMNS, I.E. THE
C PRODUCT REPRESENTATION IN TRANSPOSITIONS OF THE
C PERMUTATION WHICH WAS APPLIED TO ROWS AND COLUMNS
C OF A IN THE FACTORIZATION PROCESS
C TRAC IS A RESULTANT ARRAY OF SUBROUTINE MFSS
C TRAC MUST BE OF DOUBLE PRECISION
C INC - INPUT VARIABLE WHICH SHOULD CONTAIN THE VALUE ZERO
C IF THE SYSTEM OF SIMULTANEOUS EQUATIONS IS KNOWN
C TO BE COMPATIBLE AND A NONZERO VALUE OTHERWISE
C RHS - VECTOR OF DIMENSION N CONTAINING THE RIGHT HAND SIDE
C ON RETURN RHS CONTAINS THE MINIMAL LENGTH SOLUTION
C RHS MUST BE OF DOUBLE PRECISION
C IER - RESULTANT ERROR PARAMETER
C IER = 0 MEANS NO ERRORS
C IER =-1 MEANS N AND/OR IRANK IS NOT POSITIVE AND/OR
C IRANK IS GREATER THAN N
C IER = 1 MEANS THE FACTORIZATION CONTAINED IN A HAS
C ZERO DIVISORS AND/OR TRAC CONTAINS
C VALUES OUTSIDE THE FEASIBLE RANGE 1 UP TO N
C
C REMARKS
C THE MINIMAL LENGTH SOLUTION IS PRODUCED IN THE STORAGE
C LOCATIONS OCCUPIED BY THE RIGHT HAND SIDE.
C SUBROUTINE DMLSS DOES TAKE CARE OF THE PERMUTATION
C WHICH WAS APPLIED TO ROWS AND COLUMNS OF A.
C OPERATION IS BYPASSED IN CASE OF A NON POSITIVE VALUE
C OF IRANK
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C
C METHOD
C LET T, U, TU BE THE COMPONENTS OF THE FACTORIZATION OF A,
C AND LET THE RIGHT HAND SIDE BE PARTITIONED INTO A FIRST
C PART X1 OF DIMENSION IRANK AND A SECOND PART X2 OF DIMENSION
C N-IRANK. THEN THE FOLLOWING OPERATIONS ARE APPLIED IN
C SEQUENCE
C (1) INTERCHANGE RIGHT HAND SIDE
C (2) X1 = X1 + U * X2
C (3) X2 =-TRANSPOSE(U) * X1
C (4) X2 = INVERSE(TU) * INVERSE(TRANSPOSE(TU)) * X2
C (5) X1 = X1 + U * X2
C (6) X1 = INVERSE(T) * INVERSE(TRANSPOSE(T)) * X1
C (7) X2 =-TRANSPOSE(U) * X1
C (8) X2 = INVERSE(TU) * INVERSE(TRANSPOSE(TU)) * X2
C (9) X1 = X1 + U * X2
C (10)X2 = TRANSPOSE(U) * X1
C (11) REINTERCHANGE CALCULATED SOLUTION
C IF THE SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS IS SPECIFIED
C TO BE COMPATIBLE THEN STEPS (2), (3), (4) AND (5) ARE
C CANCELLED.
C IF THE COEFFICIENT MATRIX HAS RANK N, THEN THE ONLY STEPS
C PERFORMED ARE (1), (6) AND (11).
C
C ..................................................................
C
SUBROUTINE DMLSS(A,N,IRANK,TRAC,INC,RHS,IER)
C
C
C DIMENSIONED DUMMY VARIABLES
DIMENSION A(1),TRAC(1),RHS(1)
DOUBLE PRECISION SUM,A,RHS,TRAC,HOLD
C
C TEST OF SPECIFIED DIMENSIONS
IDEF=N-IRANK
IF(N)33,33,1
1 IF(IRANK)33,33,2
2 IF(IDEF)33,3,3
C
C CALCULATE AUXILIARY VALUES
3 ITE=IRANK*(IRANK+1)/2
IX2=IRANK+1
NP1=N+1
IER=0
C
C INTERCHANGE RIGHT HAND SIDE
JJ=1
II=1
4 DO 6 I=1,N
J=TRAC(II)
IF(J)31,31,5
5 HOLD=RHS(II)
RHS(II)=RHS(J)
RHS(J)=HOLD
6 II=II+JJ
IF(JJ)32,7,7
C
C PERFORM STEP 2 IF NECESSARY
7 ISW=1
IF(INC*IDEF)8,28,8
C
C CALCULATE X1 = X1 + U * X2
8 ISTA=ITE
DO 10 I=1,IRANK
ISTA=ISTA+1
JJ=ISTA
SUM=0.D0
DO 9 J=IX2,N
SUM=SUM+A(JJ)*RHS(J)
9 JJ=JJ+J
10 RHS(I)=RHS(I)+SUM
GOTO(11,28,11),ISW
C
C CALCULATE X2 = TRANSPOSE(U) * X1
11 ISTA=ITE
DO 15 I=IX2,N
JJ=ISTA
SUM=0.D0
DO 12 J=1,IRANK
JJ=JJ+1
12 SUM=SUM+A(JJ)*RHS(J)
GOTO(13,13,14),ISW
13 SUM=-SUM
14 RHS(I)=SUM
15 ISTA=ISTA+I
GOTO(16,29,30),ISW
C
C INITIALIZE STEP (4) OR STEP (8)
16 ISTA=IX2
IEND=N
JJ=ITE+ISTA
C
C DIVISION OF X1 BY TRANSPOSE OF TRIANGULAR MATRIX
17 SUM=0.D0
DO 20 I=ISTA,IEND
IF(A(JJ))18,31,18
18 RHS(I)=(RHS(I)-SUM)/A(JJ)
IF(I-IEND)19,21,21
19 JJ=JJ+ISTA
SUM=0.D0
DO 20 J=ISTA,I
SUM=SUM+A(JJ)*RHS(J)
20 JJ=JJ+1
C
C DIVISION OF X1 BY TRIANGULAR MATRIX
21 SUM=0.D0
II=IEND
DO 24 I=ISTA,IEND
RHS(II)=(RHS(II)-SUM)/A(JJ)
IF(II-ISTA)25,25,22
22 KK=JJ-1
SUM=0.D0
DO 23 J=II,IEND
SUM=SUM+A(KK)*RHS(J)
23 KK=KK+J
JJ=JJ-II
24 II=II-1
25 IF(IDEF)26,30,26
26 GOTO(27,11,8),ISW
C
C PERFORM STEP (5)
27 ISW=2
GOTO 8
C
C PERFORM STEP (6)
28 ISTA=1
IEND=IRANK
JJ=1
ISW=2
GOTO 17
C
C PERFORM STEP (8)
29 ISW=3
GOTO 16
C
C REINTERCHANGE CALCULATED SOLUTION
30 II=N
JJ=-1
GOTO 4
C
C ERROR RETURN IN CASE OF ZERO DIVISOR
31 IER=1
32 RETURN
C
C ERROR RETURN IN CASE OF ILLEGAL DIMENSION
33 IER=-1
RETURN
END